library("ggplot2")

make_pval_plot <- function(..., title, subtitle, ylab) {
  ggplot(...) +
    geom_point(alpha = 0.85, size = 0.5) + 
    expand_limits(x = 0, y = 0) +
    scale_x_continuous(expand = c(0, 0), limits = c(0, 1), breaks = 1:10 / 10) + 
    scale_y_continuous(expand = c(0, 0), limits = c(0, 1), breaks = 1:10 / 10) +
    coord_fixed() +
    theme_bw() +
    theme(legend.position = "none",
          text = element_text(size = 10), 
          plot.title = element_text(hjust = 0.5, size = 9.5), 
          plot.subtitle = element_text(hjust = 0.5, size = 9.5), 
          plot.caption = element_text(size = 8)) +
    geom_abline(slope = 1, intercept = 0, colour = "red", linewidth = 0.25) +
    geom_hline(yintercept = 0.05, linetype = "dashed") +
    geom_vline(xintercept = 0.05, linetype = "dashed") +
    labs(title = title, 
         subtitle = subtitle, 
         x = "Original P-Value", 
         y = ylab)
}

save_pval_plot <- function(path, plot) {
  ggsave(path, plot = plot, height = 4.1, width = 4.1, units = "in")
}


reanalysis_res <- read.csv("data/reanalysis-results.csv")

rd_papers <- read.csv("data/cl-rd-papers-meta.csv")
col_to_merge <- c("estimate_code", "bw_select_strat", "p_val")
reanalysis_res <- merge(
  reanalysis_res,
  subset(rd_papers, data_avail, col_to_merge),
  by = "estimate_code"
)
rm(rd_papers, col_to_merge)


# Figure S1

figS1 <- make_pval_plot(
  subset(reanalysis_res, bw_select_strat != "cct mse optimal"),
  aes(x = p_val, y = cct_robust_pval),
  title = "Original versus CCT Robust P-Value",
  subtitle = "(Non-CCT Studies with Replication Data)",
  ylab = "CCT Robust P-Value"
)
save_pval_plot("output/figures/figureS1-pval-comp.pdf", figS1)


# Figure S2

figS2 <- make_pval_plot(
  reanalysis_res,
  aes(x = p_val, y = cct_robust_pval),
  title = "Original versus CCT Robust P-Value",
  subtitle = "(All Studies with Replication Data)",
  ylab = "CCT Robust P-Value"
)
save_pval_plot("output/figures/figureS2-pval-comp-all.pdf", figS2)


# Figure S9: RDHonest p-value comparison plot

figS9 <- make_pval_plot(
  subset(reanalysis_res, !is.na(rdhonest_tstat)),
  aes(x = p_val, y = rdhonest_pval),
  title = "Original versus RDHonest P-Value",
  subtitle = "(Studies with Replication Data Reanalyzed with RDHonest)", 
  ylab = "RDHonest P-Value"
)
save_pval_plot("output/figures/figureS9-rdhonest-pval-comp.pdf", figS9)
